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Abstract 



The use of energy functionals based on density as the basic variable is 
advocated for ab initio molecular dynamics. It is demonstrated that the 



^ . constraint of positivity of density can be incorporated easily by using square 

o, 

^r , root density for minimization of the energy functional. An ad hoc prescription 

t — 

f— ^ . for including nonlocal pseudopotentials for plane wave based calculations is 

On ■ proposed and is shown to yield improved results. Applications are reported 



for equilibrium geometries of small finite systems, viz. dimers and trimers of 



1-^ ■ simple metal atoms like Na and Mg, which represent a rather stringent test 

a: 

O ■ for approximate kinetic energy functionals involved in such calculations. 



PACS Numbers : 71.10, 31.20G, 02.70N, 36.40 



I. INTRODUCTION 

First principle density functional based molecular dynamics (DF-MD), initiated by Car 
and Parrinello (CP) [|l],3, has become a powerful technique for ab initio investigations of 
a number of properties of clusters and extended systems 0. This technique, which unifies 
the conventional density functional theory (DFT) with classical molecular dynamics, views 
the problem as that of total energy minimization involving electronic and ionic degrees 
of freedom. This is achieved via simulated annealing implemented through Langrangian 
equations of motion, which are fictituous for electronic degrees of freedom and real for ionic 
coordinates. The method has also given impetus to the development of better and faster 
techniques for large scale electronic structure calculations with fixed geometry, e.g. clusters 
involving large number of atoms [^ . The applications based on the CP formalism fall into 
two broad categories: 

1. ab initio prediction of ground state properties like equilibrium geometry, and 

2. finite temperature properties obtained via trajectories of the system moving on Born- 
Oppenheimer (BO) surface. 

A majority of the applications belong to the first category. Inspite of a number of techni- 
cal advances like accelerated algorithms involving real space approach |^, preconditioned 
conjugate gradient minimization methods P], and analytically continued energy functional 
0, it is clear that the method becomes prohibitively expensive for large scale calculations. 
Typically these algorithms scale as 0{N^Nf,) where No is the number of orbitals and Nf, is 
the number of basis functions, N^Nf, being the dominant cost of orthogonalization of elec- 
tronic orbitals. In addition, for simple metal systems, a tricky problem of charge sloshing 
0] has also been noticed which limits the timestep that can be used in dynamics. Recently, 
there have been a few attempts towards obtaining linear scaling, e.g. a non-orthogonal lo- 
calized basis formulation , finite difference real space discretization coupled with recursion 
method [H and density matrix based formulation. 



In the present work, following the Hohenberg-Kohn theorem [|TO|,|TT|, we advocate a rather 
simple method in which the total energy functional is written in terms of density as the basic 
variable. Such an orbital free method (OFM) has been used by Pearson et al [|T^. They have 
used this method for calculating equilibrium lattice separation, bulk modulus and vacancy 
formation energy of solid Sodium. The method has also been applied to the calculation of 
free energies of vacancies [|l^. As pointed out in their work, at least for simple metallic 
systems the method scales linearly with the number of atoms Na and is capable of treating 
large simulation timesteps. The accuracy of this method hinges upon the correct description 
of the kinetic energy functional. Usually the kinetic energy functional is taken to be the 
Thomas- Fermi type with the appropriately scaled Weizsacker correction. This can be further 



improved by taking into account linear response properties, leading to the Perrot form [[13 
However, there are a few points which must be critically examined. Firstly, the usual DF- 
MD methods are implemented using first principles pseudopotentials which are necessarily 
nonlocal. So far, the available orbital free formulation [|13,|13 uses local potentials only. 
Secondly, during the minimization process the positivity of the charge density p(r) must be 
strictly maintained. Thirdly, the accuracy of the results, at least so far as the bondlengths 
are concerned, which are crucial for ground state geometries and other structural properties, 
should be thoroughly assessed. 

The present work is motivated by a desire to investigate these questions. In this work, 

1. We propose an ad hoc prescription for incorporating the nonlocal contribution of the 
pseudopotential, which is missing from the earlier work. 



2. We use square root density J p{r) as the basic variable to incorporate the positivity 
constraint on density. 

3. We present results of applications of this method to the equilibrium bondlengths of 
dimer and trimer systems. This would be a stringent test as compared to the applica- 
tions to extended systems. 



We believe that a combination of density based orbital free MD and Kohn-Sham (KS) orbital 
based MD may yield a cost effective way for performing geometry optimization for large 
clusters. This can be achieved by first obtaining approximate low temperature structures 
by the present method which is 0{Na) and then performing full MD or a fast quench. As 
the kinetic energy functional improves, this way of geometry optimization should turn out 
to be a computationally tractable alternative for large clusters consisting of more than few 
hundred atoms. It may also be used for investigating the thermodynamic properties of large 
scale systems where the conventional methods may be prohibitive in terms of computer 
time. The method will be most useful for simple metal systems provided it yields acceptable 
bondlengths. It is hoped that this kind of work will give impetus to formulating better 
kinetic energy functionals. 

In the next section we give our formalism and computational details. In section III we 
present the results for Na2, Mg2 and Mgs and compare them with full self consistent KS 
calculations. 

II. FORMALISM AND COMPUTATIONAL DETAILS 

A. Total Energy Calculation 

The total energy of a system of A'^e interacting electrons and Na atoms, according to the 
Hohenberg-Kohn theorem, can be uniquely expressed as a functional of the electron density 
p(r) under an external field due to the nuclear charges at coordinates R„. 



E[p, {R„}J = T[p] + EM + Ec[p] + E,^t[p, {R„}J + £;,,({R„}j (1) 

where 

is the kinetic energy functional. The first term in this functional is the Thomas-Fermi 
term, exact in the limit of homogenous density, and the second is the gradient correction 



due to Weizsacker. It has been pointed out that instead of A = 1 which is the original 
Weizsacker value, A = | and other empirical values turn out to yield better results [|I4[. The 
exchange- correlation energy in the usual local density approximation is given by 

E,,[p\ = j p{v)e^,{p)d'r (3) 

where exc[p\ denotes the exchange and correlation energy per particle of a uniform electron 
gas of density p. 

is the electron-electron Coulomb interaction and 



E. 



ext 



p,{R„} = V{v)p{v)dh (5) 



is the electron-ion interaction where V{t) is the external potential and the last term in the 
Eq. (1), Eii, denotes the ion-ion interaction energy. 

The external potential is usually taken to be a convenient pseudopotential, and in general, 
it is nonlocal. Let us recall that in the standard pseudopotential formulation |]TB[ the nonlocal 
contribution to the electron-ion energy for an electron in the state \E'(r) is given by 

Enim = J ^{r)VpsA^)m{r)£r (6) 

where Vps^r) is the /-dependent part of the pseudopotential and Pi is the angular momentum 
projection operator. In analogy, we suggest the following expression for nonlocal contribution 
to the total energy in terms of square root density. Let 

Eni[p] = J p{r)VpsA^)PiPi-^)d'r (7) 

where 



P(r) = VP(r)- (8) 

The exchange-correlation energy is calculated using the Ceperley-Alder exchange potential 



as interpolated by Perdew-Zunger |16|. Although in the present work we use a simple 



kinetic energy (KE) functional, improved KE functionals useful at least for simple metals 
have been reported. For example, Smargiassi and Madden [Q have investigated a family 
of KE functionals, with applications to Na and Al, giving comparable accuracy as obtained 
by KS method. 

B. Dynamics 

Typically the MD procedure proceeds via two steps. The first step is to obtain the ground 
state energy E and density p(r) for a fixed geometry, which could be done via CP dynamics 
performed on electronic degrees of freedom only. However, gradient based minimization 
techniques have also been found to be effective in locating the minimum of general functions 



and are known to have fast convergence [T^. We have applied the conjugate gradient 



(CG) algorithm [p!^-|20[] for minimizing the total energy functional for a fixed geometry 
configuration. Thus, starting with a trial r-space charge density p^, the CG algorithm 
proceeds as follows: the new charge density for the k + 1*^ iteration is constructed by 
linearly combining the charge density p'^(r) and direction d'^{r) for the A;*'' iteration as 

p''+\r) = p\r) + a''d''{r) (9) 

The initial search direction is taken to be the steepest descent direction 

d\r) = -g'{r), (10) 

and successive conjugate directions (i'^(r) are defined as 

d'{r) = -g\r) + p'd'-\r), (11) 

where 

and 



P' = -4l^i- (13) 
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The search parameter a^ > is chosen to minimize the functional E\a) = E[p^ + adi'] for 
a given p^ and d!' . In the present case, we perform this hne minimization numerically, using 
Brent's method |TS], and we notice that it does not turn out to be very expensive (about 7 



function evaluations to locate an a^). It is well known that subsequent minimizations along 
the CG directions tend to introduce errors in the calculation due to finite precision. We have 
found it advantageous to restart the CG search after every few iterations with a steepest 
descent direction. To incorporate the positivity constraint on density into the minimization 
procedure, we vary p(r) (Eq. (8)) rather than p(r), with Eq. (9)-(12) reexpressed in terms 
of p(r). 

After performing minimization to a desired degree of convergence, trajectories of ions 
and fictitious electron dynamics are simulated using Langrage's equations of motion with 
Verlet algorithm. To simulate the motion on the BO surface, we start with the Lagrangian 
defined by 

L = K, + Ka-E\p,{Yin]\ (14) 



where 



Ke = p 'p{v)d?r (15) 



and 

Ka = \Y.Mn\Rn{t)\^ (16) 

is the kinetic energy of electrons and the kinetic energy of the ions placed at {R„} respec- 
tively. The dot denotes time derivative and p and M„ are respectively the fictitious mass 
of electrons and the mass of the n-th atom. The fictitious mass of the electrons /i is a 
parameter to model the classical motion of density analogous to atomic motion. The above 
Lagrangian leads to the following equations of motion 



Mjin{t) = -VnE (18) 

for electrons and ions respectively, subject to the constraint 

Ne = /p2(r)dV. (19) 



The dynamics being conservative the grand total energy Egt is a conserved quantity of 
motion. Thus, 

Egt = E, + Ea + E [p, {R„}] (20) 

is the sum of the fictitious kinetic energy of electrons, kinetic energy of ions and total 
energy of electrons. The grand total energy is the parameter used to monitor and judge 
the quality of DF-MD numerical simulations. The greatest advantage of using density as 
the variational parameter is that the constraint of orthonogonality of wavefunctions can be 
done away with. This saves considerable amount of computation, thus making the dynamics 
fast. An alternative to CP algorithm for dynamical simulations of ionic systems has been 
suggested by Payne et al [^ . In this method the electronic degrees of freedom are relaxed 
to the instantaneous ground state at the new ionic coordinates. We have used this CG- 
MD procedure along with the density predictor method [^ for dynamical simulations. The 
density in successive timesteps is constructed using the first order density predictor as 

p/({R„(t + 1)}) = p({R„(t)}) + [p({R„(t)}) - p({R„(t - 1)})] (21) 

where pi denotes a trial density which is then used for further minimization. This reduces 
the number of CG minimization steps by a factor of 2 (typically 2 or 3 CG steps have been 
found to be sufficient for convergence). 

The calculations for Na2, Mg2 and Mgs have been performed on a periodically repeated 
unit cell of length 35 a.u. with a 48 x 48 x 48 mesh. The square root charge density is 
expanded in terms of plane waves as 

8 



p(r)=^p(G)e^«-. (22) 

G 

CG dynamics involves the calculation of the first derivatives of energies with respect to p(r). 
The electrostatic energy, gradient correction to the kinetic energy and the nonlocal energy 
and their respective derivatives were calculated in Fourier space and then transformed to r- 
space. The MD calculations were performed in the conventional CP technique and conjugate 
gradient CP technique. Local calculations were conveniently performed in the r-space. The 
energy cutoff used for local calculations was 30 Rydbergs and that for nonlocal energy 
calculations 10 Rydbergs. The dimer dynamics with conjugate gradient CP technique is 
quite stable with one CG step after adjusting the density by the predictor method for 
timestep of 10 a.u. and // = 600 a.u. Calculations with a timestep of 50 a.u. required 4 CG 
steps for the trajectories to remain on the BO surface. 

III. RESULTS AND DISCUSSION 

In this section, we present our results for the bondlengths of Na2, Mg2, Mg3 and compare 
them with that obtained via conventional ab initio MD. We have performed both fixed 
geometry minimization as well as simulated annealing MD. All the results presented here 
are obtained with energy convergence upto 10^^^ for total energy minimization 

Table I shows the comparison of bondlengths via the OFM with local and nonlocal pseu- 
dopotentials and KS self consistent formulation. It is gratifying to note that the maximum 
error in the bondlength with nonlocal pseudopotential is of the order of 2%. It can be 
seen that nonlocality improves the bondlengths considerably. However, this is at the cost of 
additional operations which go as Na x Nf,. 

In order to understand the total energy behaviour using OFM and also the effect of 
nonlocality, we have plotted the total energy as a function of interatomic separartion in 
Fig. 1. The curve labeled A corresponds to the KS local results, B corresponds to the KS 
nonlocal, C and D represent the corresponding results obtained by the present method. The 
energy and the distances are measured with respect to the equilibrium quantities obtained 

9 



by respective calculations. It can be seen that although the magnitude of changes is rather 
small, the effect of nonlocality is somewhat drastic in OFM. This leads to a faster approach 
to equilibrium with the inclusion of nonlocality. However, it also shows a worrisome feature, 
namely, the effect of nonlocality in the OFM is considerably more as compared to the full 
KS calculation. To assess the quality of the charge density obtained via the OFM, we have 
plotted it in Fig. 2 along with that due to the KS method. As can be seen, there is an overall 
agreement in the nature of the density curve. The OF density is overestimated in the region 
away from the atomic site and underestimated at the atomic site. This is understandably due 
to the inexact formulation of kinetic energy functionals. However, the long range behaviour 
is identical in both the cases. The variation of total energy E (continuous line) and Eqt 
(dashed line) with time during the free MD simulation run for Mg2 is shown in Fig. 3. It 
can be seen that the grand total energy Egt is constant throughout the simulation to an 
order of 10~^. These results have been obtained by CG-MD using density predictor method. 
We have observed that for the dimer case a much higher timestep ~ 50a. u. also keeps the 
system on Born-Oppenheimer surface with a slight increase in the computational cost. It 
may be mentioned that the Langragian dynamics (CP Technique) is much faster and our 
results for a timestep ~ 20 are more or less identical with the CG dynamics results. 

IV. CONCLUSION 

In this work we have presented a fast but approximate density based ab initio MD and 
demonstrated that the bondlengths for dimers and trimers are obtained to within an accuracy 
of 2% for nonlocal and ~ 10% for local calculations. We have shown that the CG technique 
in conjugation with a simple density predictor method allows us to use large timesteps for 
dynamics. The fixed geometry minimization using CG can be obtained to a high degree 
of accuracy within less than 200 iterations. Yet another alternative to incorporate the 
positivity constraint is by constraining variation of density during the minimization (Eq. 
(9)) by restricting a'' to appropriate range |]2D|. Our preliminary investigation indicates 
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that this technique is much faster compared to the square root density minimization. This 
is understandable, because it is well known that CG works best for quadratic and near- 
quadratic functionals, and the degree of nonlinearity of the energy functional is reduced if 
expressed in terms of p rather than p. 

We believe the present results on small systems to be a stringent test and the method 
to be a viable alternative for calculating finite temperature properties of systems involving 
a large number of atoms. 
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TABLE I. Equilibrium bonglength obtained via OEM local formulation (column 2) and 
OEM nonlocal formualtion (column 3) compared with that of KS self consistent method 
(column 4) (all values in atomic units). 



System 


OEM local 


OEM nonlocal 


KS nonlocal 


Na2 


5.34 


5.62 


5.68 


Mg2 


5.35 


6.32 


6.26 


Mg3 


5.44 


5.88 


5.82 
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FIGURES 



FIG. 1. The behaviour of total energy near equihbrium separation as a function of 
distance. All distances and energies are measured with respect to their corresponding equi- 
librium values. Curve A: KS local, B: KS nonlocal, C: OFM local, D: OFM nonlocal. 



FIG. 2. charge density r^p(r) along the symmetry axis for Mg2. The dashed line corre- 
sponds to the KS result and the continuous line represents the OFM result. 



FIG. 3. Time evolution of the total energy E (solid line) and grand total energy Egt 
(dashed line) of Mg2 during a 1000 step free OFM MD run with timestep of 10 a.u. 
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